Random variables#
What you need to know
Summing independent random variables results in another random variable called sumple sum. The mean of the sample sum is different from the population mean or expectation which is an exact quantity we want to approximate by sampling.
The Law of Large Numbers is a principle that states that as the number \(N\), the sample mean approaches the population mean with a standard deviation falling off as \(N^{-1/2}\)
The Central Limit Theorem (CLT) tells us that summing independent and identically distributed random variables with well-defined means and variances results in Gaussian distribution regardless of the nature of a random variable.
A model of random walk describes the erratic, unpredictable motion of atoms and molecules, providing a fundamental model for diffusion processes and molecular motion in fluids. The number of steps to the right (or left) of a 1D random walker results in a binomial probability distribution. Following CLT binomial distribution in the large N limit can be shown to be well approximated by gaussian with the same mean and variance.
Introducing random variables#
A random variable X is a variable whose value depends on the realization of experiment or simulations.
\(X(\omega)\) is a function from possible outcomes of a sample space \(\omega \in \Omega\).
For a coin toss \(\Omega={H,T}\) \(X(H)=+1\) and \(X(T)=-1\). Every time the experiment is done, X returns either +1 or -1. We could also make functions of random variables, e.g., every time X=+1, we ear 25 cents, etc.
Random variables are classified into two main types: discrete and continuous.
Discrete Random Variable: It assumes a number of distinct values. Discrete random variables are used to model scenarios where outcomes can be counted, such as the number of particles emitted by a radioactive source in a given time interval or the number of photons hitting a detector in a certain period.
Continuous Random Variable: It can take any value within a continuous range. These variables describe quantities that can vary smoothly, such as the position of a particle in space, the velocity of a molecule in a gas, or the energy levels of an atom.
Random numbers in python#
The numpy.random has the fastest random number generators based on low-level code written in C.
The Scipy.stats has an extensive library of statistical distributions and tools for statistical analysis.
First, we take a look at the most widely used random numbers of numpy, also called standard random numbers. These are rand (uniform random number on interval 0,1) and randn (stnadard average random number with 0 mean and 1 variance).
When running code that uses random numbers results will always differ for every run. If you want code to reproduce the same result, you can fix the seed to get reproducible results:
np.random.seed(8376743)To convert random variables to probability distributions we need to generate large enough sample then perform histogramming via
np.histor directly histogram and visualize by one shot viaplt.hist()
import numpy as np
import matplotlib.pyplot as plt
X = np.random.rand(50)
print(X)
plt.plot(X, '-o')
[0.01649061 0.94809928 0.53167608 0.70297179 0.24938418 0.644101
0.63191325 0.26702288 0.54860171 0.9113438 0.42463562 0.64045427
0.5836673 0.87348577 0.87839506 0.30074767 0.99802114 0.8662603
0.59231173 0.92788546 0.48898307 0.46747427 0.90064271 0.22878622
0.50978301 0.18799548 0.80011062 0.19234139 0.37394562 0.75512087
0.90381803 0.57959818 0.62468977 0.69036726 0.73657662 0.31476815
0.51585703 0.55235361 0.03559257 0.95396299 0.75260844 0.9636808
0.92556659 0.38376044 0.94806819 0.51243658 0.5512443 0.00192689
0.69598612 0.21691083]
[<matplotlib.lines.Line2D at 0x7feac5bb1040>]
Probability Distribution of a Random Variable#
For any random variable \( X \), we are interested in finding the probability distribution over its possible values \( x \), denoted as \( p_X(x) \).
It is important to distinguish between:
\( x \), which represents a specific value the variable can take (e.g., \( 1,2, \dots, 6 \) for a die).
\( X \), which is the random variable itself, generating values \(x\) according to the probability distribution \(p(x)\).
Histogramming#
Histograms provide an empirical estimate of distributions.
Continuous distributions require density functions (PDFs), while discrete distributions use probability mass functions (PMFs).
The bin width in histograms affects visualization, especially for discrete data.
counts, edges = np.histogram(X, range=(0,1), bins=20)
print(counts, edges)
plt.hist(X, density=True)
[3 0 0 2 3 1 2 2 1 2 5 5 4 2 2 2 1 3 7 3] [0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65
0.7 0.75 0.8 0.85 0.9 0.95 1. ]
(array([0.60235264, 0.40156842, 1.00392106, 0.60235264, 0.60235264,
2.00784212, 1.20470527, 0.80313685, 0.80313685, 2.00784212]),
array([0.00192689, 0.10153631, 0.20114574, 0.30075516, 0.40036459,
0.49997401, 0.59958344, 0.69919287, 0.79880229, 0.89841172,
0.99802114]),
<BarContainer object of 10 artists>)
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson
# Generate data for continuous distribution (Normal)
np.random.seed(42)
x_continuous = np.random.normal(loc=0, scale=1, size=1000)
# Generate data for discrete distribution (Poisson)
x_discrete = np.random.poisson(lam=3, size=1000)
# Define x values for theoretical curves
x_cont_range = np.linspace(-4, 4, 1000)
x_disc_range = np.arange(0, 10)
# Plot histograms
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Continuous distribution (Normal)
axes[0].hist(x_continuous, bins=30, density=True, alpha=0.6, color='b', edgecolor='black', label="Histogram")
axes[0].plot(x_cont_range, norm.pdf(x_cont_range, loc=0, scale=1), 'r-', lw=2, label="PDF")
axes[0].set_title("Continuous Distribution (Normal)")
axes[0].set_xlabel("Value")
axes[0].set_ylabel("Density")
axes[0].legend()
# Discrete distribution (Poisson)
axes[1].hist(x_discrete, bins=np.arange(11)-0.5, density=True, alpha=0.6, color='g', edgecolor='black', label="Histogram")
axes[1].scatter(x_disc_range, poisson.pmf(x_disc_range, mu=3), color='r', label="PMF", zorder=3)
axes[1].set_title("Discrete Distribution (Poisson)")
axes[1].set_xlabel("Value")
axes[1].set_ylabel("Probability")
axes[1].legend()
plt.tight_layout()
plt.show()
Expectation and Variance#
The expectation of a random variable, \( E[x] \), represents the theoretical mean, distinguishing it from the sample mean computed in simulations.
For example, consider the difference between:
The average height of people computed from a sample of cities.
The true mean height of the entire world population.
As the sample size increases, the sample mean converges to the expectation.
Expectation can be applied to:
The variable itself (mean).
Any function of the variable (e.g., squared deviation for variance).
Expectation of a Random Variable
When \( f(x) = x \), we obtain the mean, denoted by \( \mu \):
Using the definition of expectation, we define variance, which quantifies the spread of \( x \).
Variance as the Expectation of Mean Fluctuations
We often use the shorthand notation for variance:
\(\sigma^2 = V[x]\), where \( \sigma \) is the standard deviation.
Binomial#
A an example of discrete distribution Binomial is defined by a Probability Mass Function (PMF)
\(E[n] = Np\)
\(V[n] = 4Np(1-p)\)
Random Variable
\(B(n, p)\) modeled by
np.random.binomial(n, p, size)
Show code cell source
r = np.random.binomial(n=10, p=0.6, size=2000)
fig, ax = plt.subplots(ncols=2)
ax[0].plot(r, color='blue', label='trajectory')
ax[1].hist(r, density=True, color='red', label = 'histogram')
ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')
ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()
Gaussian#
A an example of continuous distribution Gaussian is defined by a Probability Distribution Function
\(E[x] = \mu\)
\(V[x] = \sigma^2\)
Random Variable
\(N(a, b)\) modeled by
np.random.normal(loc,scale, size=(N, M))\(N(0, 1)\) modeled by
np.random.randn(N, M, P, ...)
Show code cell source
# For a standard normal with sigma=1, mu=0
r = np.random.randn(200)
fig, ax = plt.subplots(ncols=2)
ax[0].plot(r, color='blue', label='trajectory')
ax[1].hist(r, density=True, color='red', label = 'histogram')
ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')
ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()
Uniform Distribution#
A simple example of a continuous distribution is the Uniform distribution, where all values within a given range are equally likely. It is defined by the Probability Density Function (PDF):
Expectation and Variance:
\(E[x] = \frac{a + b}{2}\)
$V[x] = \frac{(b - a)^2}{12}$
Random Variable
\(U(a, b)\) is modeled by:
np.random.uniform(low, high, size=(N, M))
\(U(0,1)\) (standard uniform) is modeled by:
np.random.rand(N, M, P, ...)
Show code cell source
# For a standard uniform
r = np.random.random(200)
fig, ax = plt.subplots(ncols=2)
ax[0].plot(r, color='blue', label='trajectory')
ax[1].hist(r, density=True, color='red', label = 'histogram')
ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')
ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()
Exact vs sampled probability distributions#
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# Define range and step size
xmin, xmax, step = -4, 4, 200
dx = (xmax - xmin) / step
# Generate x values and compute normal distribution
x = np.linspace(xmin, xmax, step)
px = norm.pdf(x) # Standard normal PDF
# Check normalization
normalization = np.sum(px * dx)
print('Normalization:', normalization)
# Generate random samples
r = np.random.randn(1000) # Increased sample size for better histogram resolution
# Create the plot
plt.figure(figsize=(8, 5))
# Histogram of sampled data
plt.hist(r, bins=30, density=True, alpha=0.6, color='blue', edgecolor='black',
label=f'Sampled: mean={r.mean():.2f}, var={r.var():.2f}')
# Plot theoretical normal distribution
plt.plot(x, px, 'k-', linewidth=2, label='Exact: mean=0, var=1')
# Formatting
plt.legend(loc="upper left", fontsize=10)
plt.ylabel(r'$p(x)$', fontsize=14)
plt.xlabel(r'$x$', fontsize=14)
plt.title("Comparison of Sampled Data with Normal Distribution", fontsize=12)
plt.grid(alpha=0.3)
plt.show()
Normalization: 0.9949421840184148
Learn about Transforming Random Variables#
When a random variable \(X \) is transformed by adding, multiplying by a constant, or applying a function \( Y = f(X) \), its probability distribution changes accordingly from \( p(x) \) to \( p(y) \).
Two commonly used transformations in statistical modeling involve generating specific distributions from standard forms:
Generating a Gaussian (Normal) distribution from a standard normal:
\[ N(\mu, \sigma^2) = \mu + \sigma \cdot N(0,1) \]Generating a Uniform distribution from a standard uniform:
\[ U(a, b) = (b - a) \cdot U(0,1) + a \]
These transformations are frequently used to construct random samples from desired distributions in simulations and statistical mechanics.
Transforming Random Variables
When transforming a random variable \( X \) to a new variable \( Y = f(X) \), the probability density functions are related by a Jacobian factor to account for how the transformation stretches or compresses the distribution:
which gives:
Examples of Simple Transformations:
Addition: \( Y = X + a \)
The probability remains unchanged except for a shift:
\[ p(y) = p(x + a) \cdot 1 \]
Multiplication: \( Y = aX \)
The distribution scales with a factor \( \frac{1}{|a|} \):
\[ p(y) = p(x) \cdot \frac{1}{|a|} \]
These transformations yield useful properties:
Shifting the Mean:
\[ E[X + a] = E[X] + a \]Scaling the Variance:
\[ V[aX] = a^2 V[X] \]
Using these properties, we can generate:
A Gaussian (Normal) distribution from a standard normal:
\[ N(\mu, \sigma^2) = \mu + \sigma \cdot N(0,1) \]A Uniform distribution from a standard uniform:
\[ U(a, b) = (b - a) \cdot U(0,1) + a \]
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
# Parameters
mu, sigma = 5, 2 # Mean and standard deviation for Gaussian
a, b = 2, 8 # Bounds for Uniform
# Generate standard distributions
std_normal = np.random.randn(10000) # N(0,1)
std_uniform = np.random.rand(10000) # U(0,1)
# Transform distributions
normal_dist = mu + sigma * std_normal # N(mu, sigma^2)
uniform_dist = a + (b - a) * std_uniform # U(a, b)
# Plot Distributions
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Normal Distribution
axes[0].hist(normal_dist, bins=40, density=True, alpha=0.6, color='b', edgecolor='black')
axes[0].set_title(f"Transformed Normal Distribution N({mu}, {sigma}²)")
axes[0].set_xlabel("x")
axes[0].set_ylabel("Density")
# Uniform Distribution
axes[1].hist(uniform_dist, bins=20, density=True, alpha=0.6, color='g', edgecolor='black')
axes[1].set_title(f"Transformed Uniform Distribution U({a}, {b})")
axes[1].set_xlabel("x")
axes[1].set_ylabel("Density")
plt.tight_layout()
plt.show()
Sum of Two Random Variables#
Consider the sum of two random variables, such as:
The sum of numbers obtained from rolling two dice.
The sum of two coin flips (e.g., heads = 1, tails = 0).
Sum of kinetic eneries of ideal gas.
The sum of random variables is itself a random variable!
We want to understand how to described the properties of summed random variables as they offer a prototype of how large systems emerge froms mall components.
Given probability distirbution of \(X_1\) and \(X_2\) how do we find probability distribution of X?
Expectation and Variance of the Sum#
Expectation is always a linear operator, which follows from the definition of expectation and the linearity of integration:
However, variance is not generally a linear operator. To see this let us write explicit formula first:
Defining the mean-subtracted variables: \(Y_i = X_i - E[X_i]\) we express variance of sum in terms of variances of component random variables
Since \( V[X_i] = E[Y_i^2] \), this simplifies to:
The cross term is called Covariance which measures the degree to which two random variables vary together:
Covariance and Correlation of Two Random Variables
If \(Cov > 0 \) or \(Cov < 0 \) we have positive/negative correlation and if \(Cov=0\) the variables are uncorrelated (but not necessarily independent)
To obtain a scale-independent measure, we define the correlation coefficient:
In the special case where \(X_1 \) and \(X_2 \) are independent, covariane is zero and we have additivity of variances!
This result is fundamental in statistical mechanics, probability theory, and the sciences, as it explains why variances add for independent random variables.
Sum of \( N \) Random Variables#
Consider a sequence of independent and identically distributed (i.i.d.) random variables, \( X_1, X_2, \ldots, X_n \).
Since they are identically distributed, each variable has a well-defined mean \( \mu \) and variance \( \sigma^2 \).
Our goal is to understand how the sum and mean of these variables depend on the sample size \( n \).
Sample Sum and Sample Mean
\( S_n \) is the sample sum, and \( M_n \) is the sample mean.
These quantities fluctuate with sample size \( n \), but we expect them to converge to their expectations for large \( n \)
Because the random variables are independent, all cross terms vanish for \(i \neq j\) between mean subtracted variables \(Y_i = X_i-E[X_i]\)
Mean and Variance of the Sum of i.i.d. Random Variables
Expectation of the Sum:
Variance of the Sum:
Law of Large Numbers#
For the sample mean the result of summatiion of i.i.d variables implies
Thus, the sample mean is an unbiased estimator of \( \mu \), and its variance decreases as \( 1/n \), meaning that the estimate becomes more stable as \( n \) increases.
Law of Large Numbers (LLN)
Implication:
The sample mean provides a reliable estimate of \(\mu\) for large \(n\).
The variance of \(M_n\) decreases as , meaning fluctuations shrink as \(1/\sqrt{n}\).
This justifies ensemble averaging in statistical mechanics, ensuring macroscopic observables (e.g., temperature, pressure) are stable and predictable.
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
# Number of trials and runs
N, runs = int(1e5), 30
# Store fractions of heads for each trial in each run
fractions = np.zeros((runs, N))
# Simulate coin tosses
for run in range(runs):
# Generate coin tosses (0 for tails, 1 for heads)
tosses = np.random.randint(2, size=N)
# Calculate cumulative sum to get the number of heads up to each trial
cum_heads = np.cumsum(tosses)
# Calculate fraction of heads up to each trial
fractions[run, :] = cum_heads / np.arange(1, N+1)
# Plotting
plt.figure(figsize=(14, 8))
# Plot all runs with low opacity
for run in range(runs):
plt.plot(fractions[run, :], color='grey', alpha=0.3)
# Highlight first run
plt.semilogx(fractions[0, :], color='blue', linewidth=2, label='Highlighted Run')
# Expected value line
plt.axhline(y=0.5, color='red', linestyle='--', label='Expected Value (0.5)')
plt.xlabel('Number of Trials')
plt.ylabel('Fraction of Heads')
plt.title('Law of Large Numbers: Fraction of Heads in Coin Tossing')
plt.legend()
<matplotlib.legend.Legend at 0x7feaa79dca90>
The Central Limit Theorem (CLT)#
Central Limit Theorem asserts that the probability distribution function or PDF of sum of random variables becomes gaussian distribution with mean \(n\mu\) and \(n\sigma^2\).
Note that CLT is based on assumption that the mean and variance, \(\mu\) and \(\sigma^2\), are finite!. Thus, CLT does not hold for certain power-law distributed random variables.
Central Limit Theorem (CLT)
Sum of any i.i.d variables (even if they are not gaussian) leads to normally distributed random variable (the sum \(s_n\))
The probability density function (PDF) of \(S_n\) is approaching gaussian:
If we subtract mean and scale the sample sum by its standard deviation we will get a standard normal distribution.
Show code cell source
from scipy.stats import norm
# Number of coin tosses in each experiment, number of experiments
N, runs = 100, 1000
# Simulate coin tosses: num_experiments rows, num_tosses_per_experiment columns
tosses = np.random.randint(2, size=(N, runs))
# Calculate means of each experiment
M = np.mean(tosses, axis=0)
z = ( M-M.mean() ) / np.std(M)
# Plotting the distribution of sample means
plt.figure()
plt.hist(z, density=True, bins=30)
plt.title('Distribution of Sample Means of Coin Tosses')
plt.xlabel('Sample Mean')
plt.ylabel('Density')
zs = np.linspace(z.min(), z.max(), 1000)
plt.plot(zs, norm.pdf(zs),'k', label='mean=0, var=1')
plt.legend()
<matplotlib.legend.Legend at 0x7feaa7fb5fd0>
Example of CLT applied to random walk problem
Applying the formulas to random walk model we get mean and variance for single step
Since steps of a random walker are independent we can compute the variance of a total displacement by multiplying mean and varaince of a single step by N
The variance of the mean \(\bar{x} = x/N\) would then be:
::
Simulating a 1D unbiased random walk#
Each random walker will be modeled by a random variable \(X_i\), assuming +1 or -1 values at every step. We will run N random walkers (rows) over n steps (columns)
We then take cumulative sum over n steps thereby summing n random variables for N walkers. This will be done via a convenient
np.cumsum()method.
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
def rw_1d(n, N):
"""
Simulates a 1D symmetric random walk.
Parameters:
n (int): Number of steps.
N (int): Number of walkers.
Returns:
np.ndarray: A (n, N) array where each column represents a walker's trajectory.
"""
# Generate random steps (-1 or +1) for all walkers
steps = np.random.choice([-1, 1], size=(n, N))
# Compute cumulative sum to get displacement
rw = np.cumsum(steps, axis=0)
# Ensure the initial position is zero
rw = np.vstack([np.zeros(N), rw]) # Adds a row of zeros at the start
return rw
# Example usage: Simulate and plot a few random walks
n_steps = 1000
n_walkers = 3
rw = rw_1d(n_steps, n_walkers)
plt.plot(rw)
plt.ylabel('X (displacement)')
plt.xlabel('n (steps)')
plt.title('1D Random Walk')
plt.show()
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# Simulate 1D random walk
def rw_1d(n_max, N):
steps = np.random.choice([-1, 1], size=(n_max, N))
return np.cumsum(steps, axis=0)
# Parameters
n_max = 1000 # Time steps
N = 1000 # Number of walkers
rw = rw_1d(n_max, N)
# Define time snapshots
time_snapshots = [10, 100, 500, 900]
# Create a multi-column subplot
fig, axes = plt.subplots(nrows=2, ncols=len(time_snapshots), figsize=(15, 6))
for i, t in enumerate(time_snapshots):
# Plot random walk trajectories
ax = axes[0, i]
ax.plot(rw[:, :50], alpha=0.3) # Show 50 trajectories for clarity
ax.axvline(x=t, color='black', lw=2) # Mark current time step
ax.set_xlabel('t')
ax.set_ylabel('X')
ax.set_title(f'Time t={t}')
# Histogram of positions at time t
ax_hist = axes[1, i]
ax_hist.hist(rw[t, :], bins=30, color='orange', density=True, alpha=0.6, label=f't={t}')
# Gaussian overlay
x = np.linspace(-100, 100, 1000)
y = stats.norm.pdf(x, 0, np.sqrt(t))
ax_hist.plot(x, y, color='black', lw=2, label='Normal')
ax_hist.set_xlim([-100, 100])
ax_hist.legend()
ax_hist.set_title(f'$\sigma/t$ = {np.var(rw[t, :])/t:.3f}')
fig.tight_layout()
plt.show()
Mean square displacement (MSD) of a random walker#
After time n number of steps (or time t) how far has random walker moved from the origin?
We quantify this by computing Mean Square Displacement (MSD). Note that the mean is computed over N number of simulated trajectories (ensemble average). Invoking central limit theorem, or simply realizing that off diagonal terms drop off we end up with the same result as in LLN.
Show code cell source
n, N = 2000, 1000
rw = rw_1d(n, N)
t = np.arange(n)
R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time
msd = np.mean(R2, axis=1) # Notice we average over N
plt.loglog(t, np.sqrt(msd), lw=3)
plt.loglog(t, np.sqrt(t), '--')
plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
def rw_2d(n, N):
"""
Simulates a 2D symmetric random walk.
Parameters:
n (int): Number of steps.
N (int): Number of trajectories.
Returns:
np.ndarray: A (n+1, N, 2) array where each trajectory is stored in the last dimension.
"""
# Define possible step directions (right, up, left, down)
steps = np.array([(1, 0), (0, 1), (-1, 0), (0, -1)])
# Generate random step indices and map to step directions
random_steps = steps[np.random.choice(4, size=(n, N))]
# Prepend an initial position at (0,0) for all walkers
rw = np.zeros((n + 1, N, 2), dtype=int)
rw[1:] = np.cumsum(random_steps, axis=0) # Compute displacement over time
return rw
# Example usage: Simulate and plot first three random walkers
n_steps = 1000
n_walkers = 100
traj = rw_2d(n_steps, n_walkers)
plt.plot(traj[:, :3, 0], traj[:, :3, 1]) # Plot first three random walkers
plt.xlabel('X (displacement)')
plt.ylabel('Y (displacement)')
plt.title('2D Random Walk')
plt.show()
Show code cell source
# Re-import necessary libraries due to execution state reset
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import HTML
# Parameters
num_steps = 200 # Number of steps in the random walk
step_size = 1 # Step size
# Generate random steps in 3D
np.random.seed(42)
steps = np.random.choice([-step_size, step_size], size=(num_steps, 3)) # Random steps in x, y, z
positions = np.cumsum(steps, axis=0) # Cumulative sum to get positions
# Create figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim([positions[:, 0].min()-1, positions[:, 0].max()+1])
ax.set_ylim([positions[:, 1].min()-1, positions[:, 1].max()+1])
ax.set_zlim([positions[:, 2].min()-1, positions[:, 2].max()+1])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Random Walk')
# Initialize the line and point
line, = ax.plot([], [], [], 'b-', lw=2)
point, = ax.plot([], [], [], 'ro')
# Update function for animation
def update(frame):
line.set_data(positions[:frame+1, 0], positions[:frame+1, 1])
line.set_3d_properties(positions[:frame+1, 2])
point.set_data([positions[frame, 0]], [positions[frame, 1]])
point.set_3d_properties([positions[frame, 2]])
return line, point
# Create animation
ani = animation.FuncAnimation(fig, update, frames=num_steps, interval=50, blit=False)
# Convert animation to HTML
html_anim = HTML(ani.to_jshtml())
plt.close(fig) # Prevent duplicate display
# Display the animation
html_anim
Animation size has reached 21025445 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.